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1.0 INTRODUCTION 


An increase in range or speed of the U.S. Navy’s surface vessels, submarines, underwater 
vehicles, and weapons can be achieved by reducing the skin friction drag of these objects. 
Microbubble drag reduction (MDR) is a unique flow control technique that employs 
injection of gas into a liquid turbulent boundary layer to form microbubbles that can 
dramatically reduce the skin friction drag. This technique, which is able to provide drag 
reductions of as much as 80%, offers great potential in Naval applications. 

McCormick and Bhattacharya (1973) reported the first microbubble drag reduction 
experiments. During the past decades, many research efforts have been devoted to 
microbubble drag reduction (Merkle and Deutsch 1992). The work conducted by 
researchers in the former Soviet Union and in the United States, primarily by the Applied 
Research Laboratory (ARL) at The Pennsylvania State University, provided the 
benchmark in microbubble drag reduction research. It has been found that there are many 
factors that influence microbubble drag reduction, including air jet flow rate, injection 
process, free stream velocity, pore size, buoyancy, and surface configuration. As 
evidenced by published papers in the open literature, most of the previous studies of 
microbubble drag reduction were conducted experimentally. Due to the complexity of the 
microbubble boundary layers, theoretical investigations have fallen behind the progress 
made by experimental studies. It is recognized that a better understanding of the 
microbubble drag reduction mechanism is critical to its optimal performance with 
minimal use of gas volume in practical applications. 

In recent years, analytical computational modeling of microbubble drag reduction has 
been attempted by several researchers (Madavan et al. 1985; Marie 1987; Kim and 
Cleaver 1995; Meng and Uhlman 1998; Xu et al. 2002; Kunz et al. 2003) to reveal the 
mechanism of the microbubble drag reduction phenomenon. A study by Legner (1984) 
proposed a simple stress model for gas bubble drag reduction and indicated that the drag 
reduction was caused by a combination of density reduction and turbulence modification. 
Meng and Uhlman (1998) suggested that bubble splitting was a plausible basic 
mechanism for reducing turbulence in a microbubble-laden turbulent boundary layer. 
These efforts made impressive progress toward the in-depth understanding of the 
mechanism from various angles. However, the available theoretical work is not sufficient 
to answer all the questions associated with microbubble drag reduction. The relative 
importance of postulated mechanisms remains unclear. One issue of specific interest is 
bubble breakup during microbubble drag reduction, which can give rise to turbulence 
attenuation. It has also been suggested that a simple density effect is the dominant source 
of drag reduction. Very recent experiments at ARL have also shown that significantly 
more drag reduction can be obtained on rough surfaces than on smooth surfaces with 
microbubble drag reduction. It is not understood why this phenomenon occurs. 

In this project, the Hemispheric Center for Environmental Technology (HCET) at Florida 
International University has conducted computational studies of the microbubble drag 
reduction mechanism. The objective of this research was to assess the roles of mixture 
density variation and surface roughness in microbubble drag reduction. 


Analytical Assessments Related to Microbubble Drag Reduction 


1 




HCET-2003-M014-001 -04 


This research was carried out in collaboration with the Applied Research Laboratory at 
The Penn State University, where substantial experimental research and theoretical 
understanding and modeling research have been conducted in this area. The principal 
work scope was carried out by HCET with ARL providing data and guidance and 
recommending models and analysis tools to support the research. 
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2.0 NUMERICAL MODEL 


During this project, computational assessment of the role of mixture density variation in 
microbubble drag reduction was carried out. To perform this assessment, a two- 
dimensional computational fluid dynamics (CFD) model of microbubble-laden flow over 
a flat plate was developed. The model consisted of Reynolds-averaged Navier-Stokes 
(RANS) transport equations and a standard k-O) turbulence model (Wilcox 1998) with a 
low Reynolds number correction. This model was designed to be applied throughout the 
boundary layer when the near wall mesh is sufficiently fine. The RANS transport 
equations are 


9 / \ 9 / \ dp d 


dt 


dxj 


3x. dxj 


^ du t + duj 2^ du, ^ 
dxj dx t 3 u dx, 


+ ^-(~ Pu'u'j) ( 2 ) 


The last term in equation (2) is called Reynolds stresses and is related to mean velocity 
gradients using Boussinesq hypothesis: 


-pu i u j =fl, 


du, | duj 
dx, dx, 


> J 


pk + p, 


du, 
dx: 


ij 


( 3 ) 


The transport equations for the turbulence kinetic energy k and the specific dissipation 
rate <»are 


f-W+^ =;T- 


r — 

V kdx Jj 


+ G k -Y k 


| w4k)=|- 


r 


d(Q 

dx 


+ G a -Y„ 


U 


( 4 ) 

( 5 ) 


where T k and T m are the effective diffusivity of k and 0) , respectively; G k is the 

generation of turbulence kinetic energy due to mean velocity gradients; Geo is the 
generation of co, and Y k and T^are the dissipation of k and 0due to turbulence. 

The low Reynolds number correction is achieved by introducing damp coefficient or* into 
the turbulent viscosity equation as shown below: 


& 



( 6 ) 


* 

a 


y 1 + Rc, /R k J 


( 7 ) 
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where ar* = 1; 


Re, = —; R k = 6; 

pCQ 


a, 


-A. 

3 ’ 


Wall boundary condition for the turbulence kinetic energy is 

^ = 0 
dn 

where n is the local coordinate normal to the wall. 


fa =0.072 


( 8 ) 


Wall boundary condition for the specific dissipation rate is discussed in detail in section 
4.0, where the effect of surface roughness on microbubble drag reduction is addressed. 

Mixture density variation due to microbubbles was modeled by introducing C0 2 gas 
species and using the species transport model. In the species transport model, the local 
mass fraction of each species, y„ is predicted by solving a convection-diffusion equation 
for the t'th species: 

jj{pY,)+V-(favY t ) = -V l+fa + S t (9) 


where 7, is the diffusion flux of species i, R, is the net rate of production of species i by 

chemical reaction (/?, = 0 in this case), and S, is the rate of creation from any sources. For 
turbulent flows such as those considered here, the diffusion flux has the following form: 


Ji = 


pD +-^- vy. 

P , ’ m Sc, 


( 10 ) 


where Sc, is the turbulent Schmidt number, p., is the turbulent viscosity, and D im is the 
diffusion coefficient for species i in the mixture. 

The C0 2 gas was introduced as species mass source in the first layer of cells along the 
porous section of the flat plate (“Wall 2” on Figure 1). Zero-gradient condition for all 
species was used on the flat plate. Mixture density was computed using the volume- 
weighted mixing law. This model does not capture the physics of the microbubbles; 
however, it allows one to ascertain whether a simple mixture density variation effect is 
the dominant source of microbubble drag reduction. The model was solved using the 
FLUENT 6 CFD solver. 


A 712 mm by 250 mm computational domain (Figure 1) with 112 x 64 grid points 
(Figure 2) was used to solve the model. The height of the computational domain was 
selected so that it would be at least 20 turbulent boundary layer thicknesses, which, for 
the present case, was about 8 = 10 mm. Wall-normal clustering of cells was used to 
resolve the boundary layer, and axial clustering was used near the ends of flat plate 
sections. Dimensions of the computational domain and boundary conditions are shown 
schematically in Figure 1. Constant velocity boundary condition is used in the domain 
inlet, symmetry boundary conditions are used in the leading part of the domain as well as 
in the far field above the flat plate, no slip boundary condition is applied at the flat plate, 
and constant pressure boundary condition is applied at the domain outlet. Values of k and 
O) at the domain inlet were specified as 1.2xl0' 5 m 2 /s 2 and 1.2x10 V 1 , respectively. The 
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same values were used as initial guesses for k and (O. To compute gas inlet volumetric 
flow rates, the depth of the domain was assumed to be 102 mm. The dimensions of the 
domain corresponded to those of the Merkle-Deutsch (1992) flat plate experimental 
configuration to facilitate comparison of the results. 
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Walll 


/ Wall 2 


Wall 3 
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pressure 

outlet 


0.25 m 


0.612 m 


Figure 1. Schematic diagram of the computational domain. 



Figure 2. The 112 x 64 computational grid. 


Seven gas input flow rates and two free stream velocities were modeled to study the 
effect of mixture density variation along with the reference cases with no bubbles. The 
input conditions for computed cases are summarized in Table 1. 
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To validate that the computational grid satisfied y + < 1 at the first cell of the wall, the y + 
at the first cell was plotted along the flat plate for all computational cases. A sample plot 
for case Q3 with free stream velocity of 10.9 m/s is shown in Figure 3. The y + values 
along the flat plate for other studied cases had similar values and are not presented here. 


Table 1. 

Input conditions of computational cases 


Case 

C0 2 Volume 
Flow Rate, m 3 /s 

Nondimensional Flow 
Rate, Q/U„A 

Free Stream 
Velocity, m/s 

Re L 

based on the total 
plate length 


HMOHi 

0 

10.9 

6.67 x 10° 

Q1-V10 


0.0025 

10.9 


Q2-V10 

0.001 

0.005 

10.9 

■■QSKlEflH 


0.002 

0.01 

10.9 

6.67 x 10° | 


0.003 

0.015 

10.9 

6.67 x 10 6 

Q5-V10 

0.004 

0.02 

10.9 

6.67x10° 

Q6-V10 

0.005 

0.025 

10.9 

6.67 x 10° 

Q7-V10 

0.006 

0.03 

10.9 

6.67 x 10° 






Q0-V4 


0 

4.2 

2.57 x 10 6 

Q1-V4 


0.0025 

4.2 

2.57 x 10° 

Q2-V4 

0.0004 

0.005 

4.2 

2.57 x 10° 

Q5-V4 

0.0015 

0.02 

4.2 

2.57 x 10° 

Q7-V4 

0.0023 

0.03 

4.2 

2.57 x 10° 


In order to validate the model further, simulated boundary layer velocity profile for the 
case without gas injection was compared to standard law-of-the-wall curves. Boundary 
layer velocity profile in the outlet plain of the computational domain is plotted in Figure 
4 in inner variables and compared to the standard curves. It is seen from the figure that 
the simulated profile is in good agreement with standard curves. 
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Figure 3. The y + in the first cell of the wall along the flat plate for 
case of Q3 and U„=10.9 m/s. 



Figure 4. Comparison of simulated boundary layer velocity 
profile with standard law of the wall curves. 
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Figure 5 compares computed drag reduction with the experimental data of Merkle and 
Deutsch (1992) and with the ensemble averaged multifield two-fluid modeling results of 
Kunz et al. (2003). In this figure, the plate drag coefficient normalized by the single¬ 
phase drag coefficient is plotted versus nondimensional gas injection flow rate. The 
figure shows not only that the results of the current model are close to those of the more 
advanced model of Kunz et al. but also that the current model correctly predicts smaller 
drag reduction for lower free stream velocity as observed in the experiments. 


1 . 0 - 


0.8 4 


O 

$ 


s 


0.6-1 


0.4-1 


0.2 A 


0.0 


—FIU-HCET (numerical U=10.9 m/s) 

—T— FIU-HCET (numerical U=4.2 m/s) 

—Kunz et al. (numerical U=10.9 m/s) 

—A— Merkle-Deutsch U=12.4 m/s (experiment) 
—O— Merkle-Deutsch U=4.2 m/s (experiment) 



-i—|—i—|—i—|—i—|—i—|—i—|—i—|—i—| 

0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 


Q/UA 

Figure 5. Comparison of the computed plate drag coefficient with 
the experimental data of Merkle and Deutsch (1992) and the 
numerical model of Kunz et al. (2003). 
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3.0 DENSITY RATIO EFFECT 


A parametric study of density ratio effect on drag reduction was carried out. Density ratio 
is defined as the ratio of density of water to that of injected gas. To perform this study, 
the density of CO 2 gas was varied in the range from 1 kg/m 3 to 200 kg/m 3 ; all other 
material properties of CO 2 were not changed. Table 2 summarizes studied input 
parameter values. 


Table 2. 

Input parameters for density ratio effect studies 


Density Ratio 


Gas Flow Rate, Q/UooA 


Free Stream Velocity, m/s 


1000 

559 

333 

250 

125 

50 

30 

10 

5 


0.005 

0.02 


10.9 


In order to investigate the effect of density ratio on microbubble drag reduction, one has 
to ensure that the volumetric concentrations of injected gas are identical for different 
density ratio values. Figure 6 shows volumetric concentration profiles for different 
density ratios at two vertical profiles above the flat plate. These profiles are designated as 
“outlet” and “profile 0.4” in Figure 1. The free stream velocity and nondimensional gas 
flow rate for all the cases shown in Figure 6 were the same and equal to 10.9 m/s and 
0.005, respectively. 

To quantify the difference between the concentration profiles, they were plotted versus 
the average concentration profile along with ±5% deviations from it as shown in Figures 
7 and 8. It is seen from Figures 7 and 8 that all the profiles at the “Outlet” and “Profile 
0.4” positions are within 5% of the average profile. This indicates that presented 
simulation cases can be used to investigate the effect of density ratio on the microbubble 
drag reduction. 

In a similar fashion, Figure 9 shows gas volumetric concentration profiles for different 
density ratios at “Outlet” position for the cases with free stream velocity of 10.9 m/s and 
nondimensional gas flow rate of 0.02. Figure 10 shows concentration profiles plotted 
versus the average concentration profile along with ± 5% deviations from average profile. 
Figure 10 shows that for nondimensional gas flow rate of 0.02 gas volumetric 
concentration profiles are within 5% of the average profile, and thus the simulation cases 
can be used for density ratio effect study. 
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Figure 6. Computed gas volume fraction profiles in the boundary 
layer at two vertical profiles above the flat plate (U=10.9 m/s, 
Q/ILA = 0.005). 
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Figure 7. Comparison of gas volume concentration profiles at 
“Profile 0.4” above the flat plate (U=10.9 m/s, Q/U A = 0.005). 
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Figure 8. Comparison of gas volume concentration profiles at 
“Outlet” above the flat plate (U=10.9 m/s, Q/U_A = 0.005). 



Figure 9. Computed gas volume fraction profiles in the boundary 
layer at “Outlet” vertical profile above the flat plate (U=10.9 m/s, 
Q/U„A = 0.02). 
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Figure 10. Comparison of gas volume concentration profiles at 
“Outlet” above the flat plate (U=10.9 m/s, Q/U^A = 0.02). 


Effect of density ratio on velocity profiles in the outlet of the computational domain is 
shown in Figure 11 for cases with nondimensional gas flow rate of 0.005. As seen from 
Figure 11, the velocity profiles for all density ratios except DR=5 are practically 
identical. The velocity profile for DR=5 has slightly lower gradient for 1 mm < y < 10 
mm than the rest of the profiles. 

Effect of density ratio on the drag coefficient is shown in Figure 12. On this figure, drag 
coefficient value for DR=1 (p gas = /Wr) is taken as that of single-phase flow with no gas 
injection. From Figure 12 it is seen that for low gas injection rate (Q/U„A = 0.005) with 
increasing density ratio, drag coefficient quickly decreases from its single-phase flow 
value and remains almost unchanging for density ratios between 5 and 1000. The drag 
coefficient for DR=5 is only about 1% higher than that for DR=1000. For high gas 
injection flow rate (Q/LLA = 0.02), a gradual decrease of the drag coefficient with 
increasing density ratio is observed, which indicates that simple mixture density variation 
effect plays one of the major roles in microbubble drag reduction phenomenon. 
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4.0 SURFACE ROUGHNESS EFFECT 


Assessments of the role of surface roughness in microbubble drag reduction were 
performed using the same numerical model. Surface roughness was included in the model 
described in Section 2 using the equivalent sand grain roughness method. In this method 
roughness height in inner variables is defined as 


k: 


f 


= max 


1 . 0 , 




pK u x \ 

A J 


( 11 ) 


Next, the value of coat the wall is specified as 


0 ) = of 
w ju 


( 12 ) 


where co + in the laminar sublayer is given by 


f 


a = min 


(O' 


w * 


rJy *F 


, where of 



K< 25 
K - 25 


(13) 


Roughness heights k s of 75 pm and 150 pm were used in the simulation cases with free 
stream velocity of 10.9 m/s and nondimensional gas injection flow rate of 0.02. The input 
parameters as well as calculated drag coefficients are summarized in Table 3. 

It is seen from Table 3 that with increasing roughness height, the HCET single-phase 
model predicts increasing drag coefficient for the cases without gas injection as well as 
for the cases with gas injection. It is known from the experiments that microbubble 
concentration reaches maximum slightly above the flat plate (Merkle and Deutsch 1992), 
while near the plate the bubble concentration is almost zero. Due to this fact, HCET 
performed a series of simulations where mixture viscosity near the wall was replaced 
with water viscosity in the thin layer above the flat plate. The input conditions and result 
of these simulations are shown in the bottom rows of Table 3. The table shows that with 
increasing height of the layer where mixture viscosity is replaced by that of water, the 
drag coefficient is decreasing; however, the drag coefficient value is still higher than that 
for smooth surface case. 


Figures 13 and 14 show the effect of roughness on velocity and gas concentration profiles 
in the outlet of the computational domain. It is seen from Figure 13 that roughness effect 
on the gas concentration is very minor and confined to the boundary layer. With 
increasing roughness, the gas concentration near the flat plate is decreasing (y < 1mm), 
while further above the plate (y > 1mm) it is increasing. Roughness effect on velocity 
profiles is also small and limited to the boundary layer as seen from Figure 14. With 
increasing roughness, velocity gradient near the flat plate is decreasing. 

To quantify the amount of drag reduction achieved by introduction of microbubbles, the 
drag coefficient for the case with bubbles was nondimensionalized by its value for the 
case without bubbles obtained for the same roughness height. The results are presented in 
Figure 15. The figure indicates that increased drag reduction is obtained with increased 
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surface roughness. The same trend was observed in recent experiments at the Applied 
Research Laboratory at Penn State (Deutsch et al. 2003). 


Table 3. 

Input parameters and results of surface roughness effect studies 


Roughness 

Height, pm 

Gas Flow Rate, 
Q/ILA 

Free Stream 
Velocity, m/s 

Viscosity 

Drag 

Coefficient 

0 

0 

10.9 

water 

0.002826272 

0 

0.02 

10.9 

mixture 

0.000571554 

75 

0 

10.9 

water 

0.005405200 

75 

0.02 

10.9 

mixture 

0.000646687 

150 

0 

10.9 

water 

0.008695437 

150 

0.02 

10.9 

mixture 

0.000800603 

150 

0.02 

10.9 

water, y<0.2 mm 
mixture y>0.2mm 

0.000794396 

150 

0.02 

10.9 

water, y<0.5 mm 
mixture y>0.5mm 

0.00076992 

150 

0.02 

10.9 


0.000734048 



C 

V 

Figure 13. Surface roughness effect on gas concentration 
profiles. 
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Figure 14. Surface roughness effect on velocity profiles. 



Figure 15. Surface roughness effect on the microbubble drag 
reduction. 
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5.0 CONCLUSIONS 


Major conclusions from the performed numerical simulation study are summarized 
below: 

• Single-phase model with bubbles introduced as species mass source was able to 
predict drag reduction consistent with experimental data and more complex two- 
fluid model results. 

• For low gas injection flow rates 

o Drag reduction was predicted even for the smallest studied density ratio. 

o Drag coefficient remained almost unchanging for density ratios between 5 
and 1000. 

• For high gas injection flow rates, gradual decrease of the drag coefficient with 
increasing density ratio was predicted, which indicates that simple mixture 
density variation effect plays a very important, if not major, role in microbubble 
drag reduction phenomenon. 

• Replacing mixture viscosity with water viscosity in the thin layer above the flat 
plate caused model to predict lower drag coefficient; however, this drag 
coefficient is still higher than that of smooth surface. 

• Single-phase model correctly predicted experimentally observed increasing drag 
reduction with increasing surface roughness. 
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Nomenclature 

Symbols 

A Area, m 2 

Cf Drag coefficient, dimensionless 

DR = — Density ratio, dimensionless 

P gas 

k s roughness height, m 

p Pressure, Pa 

Q Flow rate, m 3 /s 

Re Reynolds number, dimensionless 

t Time, s 
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U, 


Ui 

u’i 


U T 


Free-stream velocity, m/s 

mean velocity components (i = 1,2, 3), m/s 

fluctuating velocity components (/ = 1, 2, 3), m/s 

friction velocity, m/s 



m + = — velocity in inner variables, dimensionless 
y wall normal coordinate, m 

+ _ P u *y wa p normal coordinate in inner variables, dimensionless 

P 

p Dynamic viscosity, Pa-s 

p Density, kg/m 3 

Tyj wall shear stress, Pa 
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